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Visual  modelling  assumes  a  graphic  form  of  model  description  and  visual  representation  of  research 
results.  The  main  element  of  modelling  system  representation  is  a  block  diagram.  Owing  to  appearance  of 
this  line  of  investigation  («object-oriented  modelling)))  a  software  was  produced  (for  example,  Simulink 
in  MatLab). 

The  examples  of  ordinary  differential  equations  solving  schemes  for  initial  stage  discharge  of  lead-acid 
battery  positive  electrode  visual  modelling  are  presented.  The  comparative  analysis  of  two  algorithms 
(a  shooting  method  and  algorithm  similar  to  back  propagation  algorithm  used  in  the  neural  networks 
theory)  also  is  carried  out. 

©  2008  Elsevier  B.V.  All  rights  reserved. 


1.  Introduction 

During  theoretical  description  of  battery  discharge  process, 
there  are  difficulties  connected  with  the  fact  that  owing  to  the 
porous  structure  the  process  proceeds  on  all  electrode  plate  thick¬ 
ness,  and  the  discharge  characteristics  casually  vary  from  point  to 
point.  It  is  also  necessary  to  take  into  account  that  in  the  lead-acid 
battery  a  non-uniform  accumulation  of  poorly  soluble  and  poorly 
conducting  reaction  product  occurs  at  the  plate  thickness,  that 
results  in  non-uniform  distribution  of  the  discharge  characteris¬ 
tics  (electrolyte  concentration,  overvoltage,  distribution  of  reaction 
product  etc.)  at  the  electrode  thickness.  As  a  rule  boundary  value 
problem  numerical  solving  is  necessary. 

Last  decades  some  new  directions  of  problem  investigation, 
including  electrochemical  problems,  have  been  appeared: 

•  visual  modelling; 

•  artificial  intelligence  methods,  including  an  artificial  neural  net¬ 
works  (ANN)  formalism  and  genetic  algorithms  [1,2].  Analogue 
of  such  biologically  motivated  approaches  is  also  new  direction, 
which  has  arisen  in  90  years  -  artificial  immune  systems  (AIS)  for¬ 
malism.  In  particular,  the  AIS  method  allows  with  good  accuracy 
to  find  local  and  global  maxima  of  multimodal  functions,  that  it 
is  practically  impossible  to  make  with  using  traditional  methods 
[3]; 
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•  application  of  calculus  mathematics  modern  systems,  including 
symbolical  algebra  systems. 

The  visual  modelling  assumes  using  graphic  form  of  model 
description  and  evident  representation  of  research  results.  In 
connection  with  this  direction  of  investigation  the  term  “object- 
oriented  modelling”  (OOM)  appeared  and  the  appropriate  software 
in  modern  calculus  mathematics  systems  was  created.  Now  the 
most  popular  visual  modelling  system  is  Simulink  (MatLab  pack¬ 
age).  A  package  VisSim  also  received  the  wide  distribution.  It  would 
be  desirable  to  note  the  developments  of  the  Russian  programmers 
(St.-Petersburg  Polytechnical  University):  AnyLogic  [4]  and  Model 
Vision  Studium  (MVS)  [5]. 

In  OOM  a  basic  element  of  simulated  system  graphic  repre¬ 
sentation  is  the  block  diagram.  The  mathematical  diagram  and,  if 
necessary,  2D  and  3D  animation  is  used  for  result  representation. 

Advantage  of  the  given  approach  is  using  the  object-oriented 
approach  including  presentation  and  modules  of  the  appropriate 
algorithms,  solving  schemes  construction  and  also  easier  under¬ 
standing  of  computing  schemes  by  the  researchers  who  are  not 
having  programmer  skills. 

The  graphic  means  of  the  appropriate  software  allow  to  observe 
visually  the  process  of  system  characteristics  calculation  and  bring 
in  the  appropriate  corrections  into  the  model.  In  electrochemical 
problems  the  researched  characteristics  could  be,  for  example,  elec¬ 
trolyte  concentration  or  electrochemical  reaction  rate  distributions 
on  the  electrode  thickness. 

In  the  article  the  examples  of  electrochemical  processes 
visual  modelling  are  considered  including  modelling  in  MATLAB 
Simulink. 
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Fig.  1.  The  scheme  of  a  Koshi  problem  (4)  and  (6)  solution  in  Simulink. 


2.  Boundary  problem  solution  in  Simulink 


Let’s  solve  problem  (2)  with  the  boundary  conditions: 


Let’s  consider  a  technique  of  drawing  up  differential  scheme 
in  Simulink  for  the  Koshi  problem  decision  arising  from  the 
description  of  such  processes  as  “reaction-diffusion”  ones.  In 
dimensionless  form  it  is  possible  to  present  governing  differential 
equation  in  a  form 

d2y/dx2  +  (p/y)  Ay /Ax  =  H2y,  ( 1 ) 

where  H  is  dimensionless  parameter  and  p  is  a  parameter  reflecting 
problem  geometry  (p  =  0, 1,  2  for  rectangular,  cylindrical,  spherical 
coordinates,  respectively). 

Atp  =  0  Eq.  (1)  has  a  form 


d2y/dx2  =  H2y, 


(2) 


where  y  is  a  dimensionless  concentration  of  liquid  reagent. 

In  particular,  Eq.  (2)  describes  quasi  steady-state  process  such 
as  “reaction-diffusion”  in  porous  media.  Dimensionless  parame¬ 
ter  H  can  be  presented  as  H  =  (L/£),  where  £  is  the  reaction  zone 
effective  depth  (i.e.  distance  from  a  surface  in  a  plate  where  the 
reaction  rate  is  appreciably  distinctive  from  zero).  Parameter  H 
depends  on  properties  of  a  plate  material,  properties  of  reagent 
solution  and  chemical  reaction  rate  inside  a  plate.  The  concept  of 
effective  reaction  zone  depth  was  introduced  by  Ya.  B.  Zeldovich  for 
the  theoretical  description  of  catalyze  process  in  a  porous  plate. 


y(0)  =  1,  ^(1)  =  0.  (3) 

It  is  supposed  that  x  =  0  corresponds  to  a  surface  and  x  =  1  —to  a 
centre  of  a  plate. 

Let’s  solve  a  problem  (2)  with  boundary  conditions  (3)  by  the 
shooting  method. 

Let’s  put  to  a  problem  (2)-(3)  in  conformity  with  a  Koshi  prob¬ 
lem: 


dy/dx  =/, 

df/dx  =  H2y,  (4) 

y(0)  =  1,  /(0)  =  a. 


where  a  is  an  unknown  parameter. 

To  find  a  parameter  a  it  is  necessary  to  solve  the  non-linear 
equation 


=  0. 


X=1 


(5) 


Herey(x,aO  is  the  values  ofy  received  as  a  result  of  the  numerical 
decision  of  system  (4)  at  fixed  value  of  a. 

Let’s  enter  notation:  dy/da  =  u ,  df/da  =  v.  Differentiating  ordinary 
differential  equation  (ODE)  system  on  parameter  a,  we  shall  receive 
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the  ODE  system: 
du/dx  =  v, 

dv/dx  =  H2u,  (6) 

u(0)  =  0,  i/(0)  =  1. 

The  scheme  of  the  Koshi  problem  decision  (4)  and  (6)  in  Simulink 
is  shown  in  Fig.  1. 

To  search  a  suitable  value  of  parameter  a  we  create  a  M-file  in 
MATLAB  with  the  text  of  the  program: 

H=10;  alpha=0.5;  fend=0;  alpha01d=l;  eps=l.*10A(-6); 
sim('Plate') 

while  (abs(alpha01d-alpha)>=eps) 

alpha01d=alpha; 

sim('Plate') 


fendCalc=xFinal(l);  vend=xFinal(2); 

alpha=alpha-(fendCalc-fend)/vend; 

end 

sim('Plate') 

Here  Plate  is  a  file  name  with  the  scheme  in  Simulink. 

For  H=  10  after  the  program  work  termination  we  receive  the 
graph  ofy(x)  shown  in  Fig.  2. 

Changing  values  of  H  we  receive  the  appropriate  decisions 
(Fig.  3)  which  are  well  agreed  with  the  given  ones  in  refs.  [6,7]  and 
also  with  the  analytical  decision  y  =  cosh[H(l  -x)]/cosh(H). 

From  Fig.  3  it  is  well  seen  that  at  small  H  the  reagent  concentra¬ 
tion  distribution  becomes  practically  homogeneous  on  all  thickness 
of  a  plate,  and  at  increase  H  the  reaction  is  replaced  in  the  superficial 
layers  of  a  plate.  In  lead-acid  batteries  the  similar  reagent  behavior 
during  the  dicharge  results  in  strongly  non-uniform  distribution  of 
a  reaction  product  on  electrode  thickness  at  a  final  discharge  stage. 

3.  Visual  modelling  lead-acid  battery  discharge  initial  stage 

As  an  example  we  consider  the  numerical  schemes  of  current 
density  and  overvoltage  initial  distribution  modelling  in  Simulink 
of  a  package  MATLAB. 


Fig.  2.  The  graph  of  a  Koshi  problem  (4)  and  (6)  solution. 

The  soluble  system  ODEs  has  a  form: 


dr]/dz  =  i/A,  (7) 

di/dz  =  g(l  -2/7).  (8) 

with  boundary  conditions 

77(0)  =  0,  (9) 

i(0)  =  0,  (10) 

i(l)  =  1,  (11) 


where  dimensionless  parameters  and  variables  are 
z  =  x/L,  i  =  i2/1 ,  77  =  rj'F/RT,  A  =  k0RT/ILF. 

There  L  is  a  half  of  electrode  plate  thickness,  m;  I  is  an  overall 
discharge  current  density,  A  m-2 ;  i2  is  a  current  density  in  a  porous 
electrode,  A  m-2 ;  rj'  is  overvoltage,  V;  R  is  an  universal  gas  constant, 
J mol-1  K;  T  is  temperature,  K;  F  is  Faraday  constant,  As  kg-eq.;  /<o 
is  the  electrolyte  conductivity  under  initial  conditions,  m-1. 

The  system  (7)— (11 )  is  a  special  case  of  more  general  ODE  sys¬ 
tem  for  the  lead-acid  battery  positive  electrode  discharge  process 
description  [8].  It  is  supposed,  that  in  the  discharge  beginning 
electrolyte  concentration,  porosity  and  electrolyte  conductivity  are 
constant,  and  the  discharge  product  is  absent. 


H  =  0.2 


Fig.  3.  The  distribution  of  the  reagent  concentration  inside  the  plates  of  various  thicknesses. 
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Fig.  4.  The  scheme  of  ODE  system  solution  for  an  initial  stage  of  the  lead-acid  battery  discharge. 


The  presence  of  three  boundary  conditions  is  caused  by  the 
fact  that  the  system  of  two  linear  equations  (1)  and  (2)  contains 
unknown  parameter  g ,  which  reflects  a  current  exchange  density 
dependence  on  phenomenological  parameters  (initial  specific  sur¬ 
face  area  of  electrode  active  material,  exchange  current  density 
etc.),  therefore  one  boundary  condition  is  used  for  numerical  value 
definition  of  parameter  g. 

The  scheme  of  the  ODE  system  decision  ( 1  )-(2)  as  Koshi  problem 
with  the  initial  conditions  (3)-(4)  is  shown  in  Fig.  4. 

With  the  help  of  given  scheme  it  is  possible  to  quickly  “play” 
various  variants  to  find  “manually”  value  of  parameter  g  close  to 
correct  one.  Having  added  two  equations  in  scheme  it  is  possible  to 
automate  parameter  selection  process.  For  this  purpose  we  add  two 
differential  equations  for  variable  u  =  dij/dg  and  v  =  di/dg  in  system 
under  consideration  by  differentiation  the  equations  (12)-(13)  on 
the  parameter g  with  the  initial  conditions  u(0)  =  0  and  i/(0)  =  1.  After 
the  decision  of  the  Koshi  problem  we  receive  the  values  of  i(  1 )  and 
y(l).  In  a  shooting  method  the  value  of  parameter  g  is  corrected 
under  the  iterative  formula  of  Newton-Rafson: 

(17) 


where  gnew  and  g0id  are  the  values  of  parameter  g  on  the  current 
and  previous  iteration  accordingly. 

In  the  simple  case  the  visual  modelling  schemes  could  be  created 
using  spreadsheets  Excel  with  the  built-in  programming  language 
Visual  Basic  for  Applications  (VBA). 

As  an  example  the  scheme  to  calculate  initial  distribution  of  cur¬ 
rent  density  and  overvoltage  in  neural  network  basis  on  a  1st  order 
Runge-Kutta  method  in  the  spreadsheets  Excel  is  partially  shown 
in  Fig.  5. 

The  neural  network  basis  represents  the  weighed  summation  of 
several  entrance  signals  operation  with  the  subsequent  transfor¬ 
mation  of  this  sum  by  means  of  formal  neuron  activation  function 
which  can  be  both  linear  and  non-linear. 

In  particular,  for  system  of  two  equations  and  the  1st  order 
Runge-Kutta  method  we  receive  the  formulas: 

yi,t+i  =  bnyht  +  buy2,t  =  0  +  ban  )yi,t  +  ha\2y2,t , 

y2,t+i  =  b2\y\,t  +  b22y2,t  =  ha21yht  +  (1  +  ha22)y2,t. 

Here  y2t,  y2t+1  is  the  decision  of  ODE  system  at  step  t  and 
t+1  accordingly.  Similar  expression  could  be  turned  out  for  a 
Runge-Kutta  method  of  higher  orders. 


Snew  —  Sold  —  [*("0  —  l]AO)> 
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We  see  that  the  expression  in  the  right  part  represents  the 
weighed  sum  with  a  weight  matrix 

^  _  f  1  +  han  hai2  \ 
y  ha2 1  1  +  ha22  )  ' 

The  scheme  is  shown  after  installation  of  problem  initial  values. 
The  scheme  enables  to  lead  search  of  parameter  on  one  of  two  algo¬ 
rithms:  shooting  algorithm  and  offered  by  us  algorithm  similar  to 
the  back  propagation  algorithm  in  ANN  formalism.  Then  minimized 
functional  will  have  a  form  (l/2)-[i(l )  —  l]2  and  we  will  receive  a 
recurrence  formula 

gnew=g0ld-y[i(l)-1Ml),  (18) 

where  y  is  a  parameter  of  training. 

The  iterations  come  to  an  end  when  the  difference  of  parameter 
g  values  on  the  current  and  previous  steps  becomes  less  than  the 
given  value. 

The  comparison  of  the  formulas  (17)  and  (18)  shows  that  it  is  to 
be  expected  the  best  convergence  of  back  propagation  algorithm 
especially  when  the  value  of  )  is  small  enough.  The  calculations 
show  that  the  algorithm  convergence  rate  is  sensitive  to  the  value 
of  parameter  y  and  in  the  calculus  error  of  value  of  i(  1)  is  less  by 
the  order  of  magnitude  in  comparison  with  the  shooting  method. 

From  Fig.  5  the  opportunity  of  parallel  calculations  is  well  seen 
during  the  numerical  decision  of  the  given  problem  when  the  sepa¬ 


rate  program  parts  could  be  carried  out  independently  on  different 
computing  devices. 

There  is  also  other  approach  to  the  decision  of  ODE  problems 
based  on  a  neurocomputing  [9]. 

4.  Conclusions 

Some  new  methods  of  the  differential  equations  systems  solving 
are  briefly  considered.  Similar  problems  are  often  occur  at  mod¬ 
elling  various  processes  including  electrochemical  ones. 

The  algorithm  of  solution  is  offered  which  is  similar  to  the  back 
propagation  algorithm  used  for  training  artificial  neural  networks. 
The  comparison  of  this  algorithm  with  shooting  one  was  made  for 
the  problem  of  lead-acid  battery  positive  electrode  discharge  mod¬ 
elling.  It  is  shown  that  the  rate  of  algorithm  convergence  depends 
efficiently  on  the  value  of  parameter  of  training  and  at  a  suitable 
choice  of  given  parameter  the  algorithm  works  faster  and  the  cal¬ 
culation  error  is  less  by  the  order  in  comparison  with  shooting 
algorithm. 
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